home *** CD-ROM | disk | FTP | other *** search
/ BCI NET 2 / BCI NET 2.iso / archives / programming / source / graphicgems4.lha / GemsIV / vec_mat / ray / Polyhedron.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-06  |  5.8 KB  |  193 lines
  1. #include "Polyhedron.h"
  2.  
  3. /************************************************************************
  4. *                                    *
  5. * Class destructor                            *
  6. *                                    *
  7. ************************************************************************/
  8.  
  9. Polyhedron::~Polyhedron()
  10. {
  11. int i;
  12.  
  13. // delete the vertices and the plane equations
  14. delete vList;
  15. delete pList;
  16.  
  17. // delete the facets
  18. for(i = 0; i < facetN; i++)
  19.     delete iList[i];
  20. delete iList;
  21. }
  22.  
  23. /************************************************************************
  24. *                                    *
  25. * This method computes the intersection between a ray and the        *
  26. * polyhedron.  It returns the distance between the origin of the ray    *
  27. * and the closest point of intersection (or 0.0 if not intersection    *
  28. * occurs).                                *
  29. * The algorithm operates as follows:                    *
  30. * For each facet of the polyhedron                    *
  31. *  + determine the intersection between the ray and the plane        *
  32. *    supporting the facet.                        *
  33. *  + if the intersection exists, project the facet on a 2D plane.    *
  34. *  + Test if the intersection point is inside the resulting polygon.    *
  35. *                                    *
  36. ************************************************************************/
  37.  
  38. double Polyhedron::intersect(vec3& ray_org, vec3& ray_dir)
  39. {
  40. int    i, j,
  41.     sNum = 0,   // Number of intersections.
  42.     axis,        // Axis along which the facet will be projected.
  43.     iNum;        // Number of intersections between an imaginary half line
  44.             // and the edges of a facet (used to determine inclusion).
  45. vec3    normal;        // Normal to the current facet.
  46. vec2    p1, p2,        // Two consecutive points on a facet.
  47.     p,        // A 2D projection of the point to test.
  48.     n;        // Normal to one of the edges of a facet.
  49. double    s[facetN],  // Space to store at most N intersections.
  50.     vd,
  51.     vn;
  52.  
  53. for (i=0; i<facetN; i++) {
  54.     // Test the case in which the ray is parallel to the facet (vd=0.0)
  55.     vd = vec3(pList[i],PD) * ray_dir;
  56.     if (vd == 0.0)
  57.     continue;
  58.  
  59.     // Find a projection axis which maximizes the area of the facet
  60.     (normal = vec3(pList[i], PD)).apply(&fabs);
  61.     if (normal[VX] >= normal[VY])
  62.     if(normal[VX] >= normal[VZ]) axis = VX;
  63.     else axis = VZ;
  64.     else if(normal[VY] >= normal[VZ]) axis = VY;
  65.     else axis = VZ;
  66.  
  67.     // Compute the point of intersection between the ray and the facet plane.
  68.     // Project it along the desired axis.
  69.     vn = pList[i] * vec4(ray_org);
  70.     s[sNum] = -vn / vd;
  71.     p = vec2(ray_org + s[sNum] * ray_dir, axis);
  72.  
  73.     // Test if the intersection point is inside the facet.
  74.     // We do this by checking for each edge of the facet if there is an
  75.     // intersection between this edge and a half-line, starting at p and
  76.     // going towards positive X. We count the number of intersections: if it
  77.     // is odd, the point is inside the facet, if not, it is outside.
  78.  
  79.     iNum = 0;
  80.     p1 = vec2(vList[iList[i][iList[i][0]]], axis);
  81.     for (j=1; j<= iList[i][0]; j++) {
  82.     p2 = vec2(vList[iList[i][j]], axis);
  83.  
  84.     // Is p in the band generated by sweeping the [p1p2] segment in the
  85.     // positive X direction ?
  86.     if ((p1[VY] - p[VY]) * (p2[VY] - p[VY]) <= 0.0) {
  87.  
  88.         // Does the half-line trivially intersect the edge ?
  89.         if (p[VX] < min(p1[VX], p2[VX]))
  90.         iNum++;
  91.  
  92.         // Does the half-line trivially miss the edge ?
  93.         else if (p[VX] < max(p1[VX], p2[VX])) {
  94.  
  95.         // tough case: compute the normal to the edge pointing towards
  96.         // positive X use the dot product between this normal and the
  97.         // p1p vector to see if the edge and the half line intersect.
  98.         // (This is similar to the Cyrus-Beck line-clipping test)
  99.         n[VX] = p1[VY] - p2[VY];
  100.         n[VY] = p2[VX] - p1[VX];
  101.         if (n[VX] < 0.0)
  102.             n = -n;
  103.         if (n * (p - p1) <= 0.0)
  104.             iNum++;
  105.         }
  106.         }
  107.     p1 = p2;
  108.     }
  109.     if (iNum % 2)
  110.     sNum++;
  111.     }
  112. return closest_intersection(s, sNum);
  113. }
  114.  
  115. /************************************************************************
  116. *                                    *
  117. * This method computes the normal vector to the polyhedron at the point *
  118. * of intersection.                            *
  119. * It does so by finding out which plane the intersection point belongs    *
  120. * to. Once this plane is determined, it simply returns the normal to    *
  121. * this plane, which is already normalized.                *
  122. *                                    *
  123. ************************************************************************/
  124.  
  125. vec3 Polyhedron::normalAt(vec3& p)
  126. {
  127. int    i,
  128.     index = 0;    // index of the facet being hit by the ray
  129. double    h,
  130.     hmin = 1e16;    // initialize hmin to some very large number
  131.  
  132. for (i=0; i<facetN; i++) {
  133.     h = fabs(pList[i] * vec4(p));
  134.     if (h < hmin) {
  135.     index = i;
  136.     hmin = h;
  137.     }
  138.     }
  139. return vec3(pList[index], PD);
  140. }
  141.  
  142. /************************************************************************
  143. *                                    *
  144. * Input from stream.                            *
  145. *                                    *
  146. ************************************************************************/
  147.  
  148. istream& operator >> (istream& s, Polyhedron& a)
  149. {
  150. int    i,j,M;
  151. mat4    T;    // local coordinates to world coordinates transformation
  152. vec3    n;    // normal to the facet
  153.  
  154. // call the implementation of the super-class
  155. s >> *((Primitive*) &a);
  156.  
  157. // create the matrix to transform local coordinates to world coordinates
  158. T = translation3D(a.pos) * (a.orient.transpose());
  159.  
  160. // read the vertices. Transform them on the fly to world coordinates
  161. s >> a.vertexN;
  162. a.vList = new vec3[a.vertexN];
  163. for (i = 0; i < a.vertexN; i++) {
  164.     s >> a.vList[i];
  165.     a.vList[i] = T * a.vList[i];
  166.     }
  167.  
  168. // read the facets.
  169. s >> a.facetN;
  170. a.iList = new int*[a.facetN];
  171. for (i = 0; i< a.facetN; i++) {
  172.     s >> M;
  173.     a.iList[i] = new int[M+1];
  174.     a.iList[i][0] = M;
  175.     for (j=1; j<=M; j++)
  176.     s >> a.iList[i][j];
  177.     }
  178.  
  179. // compute the normal and plane equation of the facet using Newell method
  180. a.pList = new vec4[a.facetN];
  181. for (i = 0; i < a.facetN; i++) {
  182.     n = vec3(0.0);
  183.     for (j = 1; j <= a.iList[i][0]; j++)
  184.     if (j != a.iList[i][0])
  185.         n += a.vList[a.iList[i][j]] ^ a.vList[a.iList[i][j+1]];
  186.     else
  187.         n += a.vList[a.iList[i][j]] ^ a.vList[a.iList[i][1]];
  188.     n.normalize();
  189.     a.pList[i] = vec4(n, - a.vList[a.iList[i][1]] * n);
  190.     }
  191. return s;
  192. }
  193.